# Make maps with states, turbines, and circles

library(sp)
library(tidyverse)
library(readxl)
library(rgeos)
library('raster')
library('geosphere')
library(maptools)
library(maps)
library(haven)
require ('rgdal')
library(Rcpp)
library('mapview') 
library(tigris)
library(usmap)
library(maps)
library(ggplot2)

# read shapefile as SpatialPolygonsDataFrame from TIGER/LINE shapefiles (tiger package)
counties <- counties(cb = FALSE, class = "sp")
counties48 <- counties[counties@data$STATEFP!="78" 
                        & counties@data$STATEFP!="66"
                        & counties@data$STATEFP!="69"
                        & counties@data$STATEFP!="72"
                        & counties@data$STATEFP!="60"
                        & counties@data$STATEFP!="15"
                        & counties@data$STATEFP!="02",]
#plot(counties48)

states <- states(cb = FALSE, class = "sp") 
#class(states)   # SpatialPolygonsDataFrame object
# check coordinate system/projection
#proj4string(states) # North American Datum 83 Geographic Coordinate System

# subset to lower 48
states48 <- states[states@data$NAME!="United States Virgin Islands"
                   & states@data$NAME!="Guam"
                   & states@data$NAME!="Commonwealth of the Northern Mariana Islands"
                   & states@data$NAME!="Hawaii"
                   & states@data$NAME!="American Samoa"
                   & states@data$NAME!="Alaska"
                   & states@data$NAME!="Puerto Rico",]
#plot(states48)

# store it as SpatialPolygons
SP.states48 <- as(states48, "SpatialPolygons")

# this data is each turbine, not just each plant: 
setwd("C:/Users/gilbe/Dropbox/EnergyImpactsAggregationBias/SteveStuff/data/lbnl")
# uswindturbines <- read_excel("uswtdb_public_070318_pid.xlsx", sheet = "Data")
uswindturbines <- read_dta("uswtdb_public_20211101_v4_2_pid.dta")

turbines48 <- subset(uswindturbines,t_state != "AK"
                     & t_state != "PR"
                     & t_state != "GU"
                     & t_state != "HI"
                     & p_tnum < 1400
                     & p_tnum > 7 # 2.5 percentile
                     & p_cap > 1 & p_year>1994)

# quantile(uswindturbines$p_tnum[uswindturbines$p_tnum<1400],probs = c(0.025))
turbcoord48 <-data.frame(cbind(turbines48$xlong, turbines48$ylat)) 
names(turbcoord48)<-c("Lon_T","Lat_T") 

# keep only relevant variables for later
turbines48<-data.frame(cbind(as.numeric(turbines48$xlong), 
                             as.numeric(turbines48$ylat),
                             turbines48$t_state, turbines48$t_county,
                             turbines48$t_fips, turbines48$p_name,
                             turbines48$p_year, turbines48$p_tnum,
                             turbines48$p_cap, turbines48$t_cap)) 
names(turbines48)<-c("Lon_T","Lat_T","t_state","County","FIPS","PName","Year",
                     "p_tnum","p_cap","t_cap") 
# convert missings to NAs
turbines48$t_cap[turbines48$t_cap==-9999] <- NA
turbines48$p_cap[turbines48$p_cap==-9999] <- NA
# could fill in average per-turbine cap, e.g., PCap/TNum if TCap = -9999
# but PCap is often missing

# Convert to SpatialPointsDataFrame
t.spoint.p <- SpatialPointsDataFrame(turbcoord48[, c("Lon_T", "Lat_T")], 
                                     data.frame(ID=seq(1:nrow(turbcoord48))), 
                                     proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

plot(states48)
points(turbcoord48,pch=20)

# transform to planar projection: 
t.spoint.mp <- spTransform(t.spoint.p,CRS( "+init=epsg:2163" )) 

c=c("ALABAMA","ARIZONA","ARKANSAS","cALIFORNIA","COLORADO", "DELAWARE",
    "DC","ILLINOIS","INDIANA","IOWA","KANSAS", "MAINE","MARYLAND","MISSOURI",
    "MONTANA", "NEBRASKA", "NEVADA","NEW MEXICO", "NORTH DAKOTA", "OKLAHOMA",
    "TENNESSEE","VIRGINIA", 'WYOMING')


namevec<-map(database="state", col="blue", fill=T,namesonly = TRUE )
map(database = "state",col = c("white", "light blue")[1+(namevec %in% tolower(c) )],fill=T)
points(turbcoord48,pch=20,cex=0.5)

# Create the circles ---------------------

# R spatial operations in meters. Keep track of number of meters in each d-mile radius
# MODIFY AS NEEDED FOR THE RADIUS OF YOUR CIRCLE - this converts meters to miles
maxdist_5 <- 8047 
maxdist_10 <- 16093.5
maxdist_20 <- 32187    
maxdist_40 <- 64374
maxdist_60 <- 96561
maxdist_80 <- 128748
maxdist_100 <- 160935

# create buffers for each each distance band
# THIS DOES NOT NEED TO BE A LOOP IF YOU HAVE ONLY ONE CIRCLE DISTANCE YOU WANT TO USE
for (k in c(20,60,100)) {
  #cbuf.k.p <- gBuffer(t.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE) # Uncomment this line if you want union of all circles, not individual circles
  assign(paste0("cbuf.",k,".p"),gBuffer(t.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE))
  assign(paste0("c.",k,".p.reproj"), spTransform(get(paste0("cbuf.",k,".p")),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
  assign(paste0("c.",k,".p.reproj"), as(get(paste0("c.",k,".p.reproj")), "SpatialPolygonsDataFrame"))  
}
cbuf.100.u <- gUnaryUnion(cbuf.100.p, id = NULL) 
plot(cbuf.100.u)
cbuf.60.u <- gUnaryUnion(cbuf.60.p, id = NULL) 
plot(cbuf.60.u)
cbuf.20.u <- gUnaryUnion(cbuf.20.p, id = NULL) 
plot(cbuf.20.u)

crs<-"+proj=longlat +datum=NAD83 +no_defs"
cbuf.100.u <-spTransform(cbuf.100.u, CRS=crs)
cbuf.100.u <- as(cbuf.100.u,"SpatialPolygonsDataFrame")

cbuf.60.u <-spTransform(cbuf.60.u, CRS=crs)
cbuf.60.u <- as(cbuf.60.u,"SpatialPolygonsDataFrame")

cbuf.20.u <-spTransform(cbuf.20.u, CRS=crs)
cbuf.20.u <- as(cbuf.20.u,"SpatialPolygonsDataFrame")

# check that it plots something
plot(states48)
plot(cbuf.100.u,col="green",add=TRUE)

# create maps to save - toggle on the "png" line and the "dev.off()" lines if you want to save
#### Figure 1 in the paper: 
#png("C:/Users/bgilbert_a/Dropbox/EnergyImpactsAggregationBias/statesturbs.png")
png("C:/Users/gilbe/Dropbox/EnergyImpactsAggregationBias/statesturbsgrey.png",width=500,height=500,units='mm',res=300)
map(database = "state",col = c("white", "light grey")[1+(namevec %in% tolower(c) )],fill=T)
points(turbcoord48,pch=20,cex=1)
dev.off()

#### First Appendix Figure: 
png("C:/Users/gilbe/Dropbox/EnergyImpactsAggregationBias/statesturbsbuffs100grey.png",width=500,height=500,units='mm',res=300)
#png("C:/Users/gilbe/Dropbox/EnergyImpactsAggregationBias/statesturbsbuffs100grey.png")
map(database = "state",col = c("white", "light grey")[1+(namevec %in% tolower(c) )],fill=T)
plot(cbuf.100.u,col="dark grey",add=TRUE)
points(turbcoord48,pch=20,cex=1)
dev.off()


### example plotting commands not used in the paper: 

#png("C:/Users/bgilbert_a/Dropbox/EnergyImpactsAggregationBias/statesturbsbuffs60.png")
map(database = "state",col = c("white", "light blue")[1+(namevec %in% tolower(c) )],fill=T)
plot(cbuf.60.u,col="green",add=TRUE)
points(turbcoord48,pch=20,cex=0.5)
#dev.off()

#png("C:/Users/bgilbert_a/Dropbox/EnergyImpactsAggregationBias/statesturbsbuffs20.png")
map(database = "state",col = c("white", "light blue")[1+(namevec %in% tolower(c) )],fill=T)
plot(cbuf.20.u,col="green",add=TRUE)
points(turbcoord48,pch=20,cex=0.5)
#dev.off()

map(database = "county",col="red")
plot(cbuf.100.u,border="black",add=TRUE)
map(database = "county",col="red")
plot(cbuf.60.u,border="black",add=TRUE)
map(database = "county",col="red")
plot(cbuf.20.u,border="black",add=TRUE)

map(database="county",regions="ILLINOIS")
plot(cbuf.20.u,border="red",add=TRUE)

map(database="county",regions="COLORADO")
plot(cbuf.20.u,border="red",add=TRUE)






# create buffers for each each distance band
# THIS DOES NOT NEED TO BE A LOOP IF YOU HAVE ONLY ONE CIRCLE DISTANCE YOU WANT TO USE
for (k in c(20,40,60,80,100)) {
  #cbuf.k.p <- gBuffer(c.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE) # Uncomment this line if you want union of all circles, not individual circles
  assign(paste0("cbuf.",k,".p"),gBuffer(c.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE))
  assign(paste0("c.",k,".p.reproj"), spTransform(get(paste0("cbuf.",k,".p")),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
  assign(paste0("c.",k,".p.reproj"), as(get(paste0("c.",k,".p.reproj")), "SpatialPolygonsDataFrame"))  
}
